Unveiling the cosmological QCD phase transition through the eLISA/NGO detector 
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We study the evolution of turbulence in the early universe at the QCD epoch using a state-of- 
the-art equation of state derived from lattice QCD simulations. Since the transition is a crossover 
we assume that temperature and velocity fluctuations were generated by some event in the previous 
history of the Universe and survive until the QCD epoch due to the extremely large Reynolds number 
of the primordial fluid. The fluid at the QCD epoch is assumed to be non- viscous, based on the fact 
that the viscosity per entropy density of the quark gluon plasma obtained from heavy-ion collision 
experiments at the RH1C and the LHC is extremely small. Our hydrodynamic simulations show 
that the velocity spectrum is very different from the Kolmogorov power law considered in studies of 
primordial turbulence that focus on first order phase transitions. This is due to the fact that there 
is no continuous injection of energy in the system and the viscosity of the fluid is negligible. Thus, 
as kinetic energy cascades from the larger to the smaller scales, a large amount of kinetic energy is 
accumulated at the smallest scales due to the lack of dissipation. We have obtained the spectrum 
of the gravitational radiation emitted by the motion of the fluid finding that, if typical velocity and 
temperature fluctuations have an amplitude (Aw)/c > 10 -2 and/or AT/T C > 10 -3 , they would be 
detected by eLISA at frequencies larger than ~ 10 -4 Hz. 
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I. INTRODUCTION 

During its evolution the universe passed through sev- 
eral phase transitions. The theory of quantum chromo- 
dynamics (QCD) predicts a phase transition in which a 
hot quark-gluon unconfined phase is converted, as the 
Universe expands and cools, into a confined hadronic 
phase. The best quantitative evidence for such transition 
is found in lattice gauge theory of QCD, which shows that 
this transition occurred at temperatures around 150—200 
MeV [H,0|- According to recent results, the phase transi- 
tion for small chemical potentials (condition expected at 
this epoch for standard cosmology and particle models) 
is merely an analytic transition (or crossover) Q. 

Since the QCD transition occurred before the 
radiation-matter decoupling, the only way to obtain 
information about it is possibly through gravitational 
waves created at this epoch. With the eLISA/NGO (New 
Gravitational wave Observatory) Q launch planned for 
the next years and new detectors being projected, e.g. 
Big Bang Observatory (BBO) and TOBA @, we expect 
gravitational radiation to be a new source of information 
about the QCD transition. 

Most previous studies on this topic have focused on 
the effects of a first order transition. In 1984, Witten 
@ considered the possibility that the QCD phase transi- 
tion may have produced a detectable gravitational wave 
signal if the transition is first order and leads to violent 
bubble collisions. Later, Applegate and Hogan Q made 
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a detailed study of the influence of the QCD transition in 
the standard cosmological model and studied the relics 
produced at that period. Polnarev [8j, demonstrated that 
polarization and anisotropy in the cosmic background ra- 
diation can be induced by gravitational waves depending 
on their characteristic length (see also @). The study 
of gravitational waves in cosmological phase transitions 
had a peak in the early 1990s. Based on the assumption 
of a first order transition, several studies were carried 
out about the nucleation, growth and collision of bub- 
bles and their relation to the generation of gravitational 
waves (see e.g. [Tol - [l2T |L A decade later, works on pri- 
mordial turbulence continued to be conducted. For ex- 
ample, Kosowsky, Mack and Kahniashvili [13[ calculated 
the stochastic background of gravitational radiation aris- 
ing from a period of cosmological turbulence, using a 
simple model of isotropic Kolmogorov turbulence pro- 
duced in a cosmological phase transition. Also, Dolgov 
and Grasso [14J showed that an inhomogeneous cosmolog- 
ical lepton number may have produced turbulence in the 
primordial plasma when neutrinos entered the (almost) 
free-streaming regime. This effect may be responsible for 
the origin of cosmic magnetic fields and give rise to a 
detectable background of gravitational waves. 

More recent works have addressed various possibilities 
for early-universe physics leading to detectable cosmolog- 
ical gravitational wave backgrounds and have analyzed in 
more detail the turbulent spectrum as well as the spec- 
trum of the gravitational wave signal fl5T - l2l| . However, 
to the best of our knowledge, there are no studies of tur- 
bulence at the QCD epoch assuming a crossover tran- 
sition and employing a realistic lattice QCD equation 
of state (EoS) for the primordial fluid. In the present 
work, we consider an EoS obtained from recent lattice 
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QCD simulations by the Wuppertal-Budapest collabo- 
ration [3j and perform hydrodynamic numerical simula- 
tions in order to obtain the turbulent spectrum as 
well as the gravitational radiation generated by the mo- 
tion of the fluid. The detectability of the inferred back- 
ground is discussed in the light of the recently published 
eLISA/NGO's sensitivity curve [11]. 

This article is organized as follows: we present firstly 
the relativistic hydrodynamic equations for a perfect fluid 
in one dimension and the here-used EoS (Sec. [TTJ) and we 
then describe the numerical method employed for the so- 
lution of these equations fSec. HIip . In Sec. lIVi we present 
the formalism to obtain the spectrum of gravitational 
waves from the hydrodynamic evolution of the fluid based 
on Ref. [23|. Finally, we present our results (Sec. IVI|) and 
conclusions (Sec. I VIII) . 



II. BASIC EQUATIONS 
A. Hydrodynamics 

To investigate the dynamics and evolution of the pri- 
mordial fluid, we should solve the equations of hydrody- 
namics in the context of general relativity. However, as 
in previous work (24j we shall assume a flat space metric 
filled with a perfect fluid whose stress-energy tensor is 
defined by 



T» v = phu^u" +pr) 



(1) 



where p is the pressure, h is the specific enthalpy defined 
as/i = l + e+ £,eis the specific internal energy and p is 
the mass density in the rest frame. 

In the one-dimensional case adopted here, the hydro- 
dynamic equations written in covariant form consist of 
three local conservation laws for the above stress-energy 
tensor and the baryon density flow J M = pu* 1 : 



V„T» V = 0, 



V M J^ = 0, 



(2) 



where we used units in which the speed of light is c = I . 
The system must be complemented with an equation of 
state having the functional form p = p(p, e). 



B. Equation of State 

In order to construct the EoS for a crossover tran- 
sition we considered the most recent results of lattice 
QCD simulations obtained by the Wuppertal-Budapest 
collaboration (see Ref. [3[), with Nf = 2 + 1, i.e. two 
light (up and down) and one heavy (strange) quarks. To 
adapt such EoS to the conditions prevailing at the QCD 
epoch in the early universe we added the contribution 
of a gas of non-interacting neutrinos (v T , v^, v e ), muons, 
electrons, photons and their antiparticles, all them de- 
scribed by an EoS of the form E' = p'/3 = gn 2 T 4 /90 
with g = 7/8 x 14 + 2 = 14.24. As showed in Ref. |[ the 
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Figure 1: EoS describing a crossover obtained using data from 
lattice QCD [3| plus the contribution of leptons and photons. 
The upper dashed line represents the Stefan-Boltzmann limit 
for a relativistic ideal gas of the same particles (psb ~ 6.77T 4 
and E S b ~ 3pss). 



addition of the charm or heavier quarks doesn't introduce 
any substantial modification in the temperature regime 
relevant for this work. Such data generate smooth curves 
for several observables which are typical of a crossover 
transition. In Fig. Q] we show the behavior of the pres- 
sure p and the energy density E — p(l + e) for a large 
range of temperatures as well as the Stefan-Boltzmann 
limit. 



III. NUMERICAL APPROACH 

The hydrodynamic equations presented in Eq. ([2]) can 
be recast as a hyperbolic system of first-order, flux- 
conservative equations of the form [25j : 



dx° ^ 



0. 



(3) 



where in the laboratory frame the state vector U con- 
tains the conserved variables (D,Si,r) written in terms 
of primitive variables w = (p, Vi,p) 



(4) 





D 




Wp 


W(w) = 


Sj 




phW 2 v 3 




T 




phW 2 -p-Wp 



and the flux vector ^F l (\v) is given by 



^(w) 



Dv l 
S 3 v l +p8) 

TV 1 + pV % 



(5) 



with W = (1 — v l v^)~ 1 / 2 being the Lorentz factor. The 
elements of the flux matrix are the mass flux Dv l , the 
momentum flux plus pressure force SjV 1 +P<5], and the 
energy flux plus pressure work tv % + pv l . Note that the 
flux vector is written in terms of the conserved variables 
of the state vector. 
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The conservation laws written in the form of a Eq. 
^ allow to solve the problem numerically through a 
wealth of methods. In the present work we shall employ a 
High- Resolution Shock- Capturing scheme (HRSC) which 
is recognized as a very efficient scheme for dealing with 
complex flows specially when discontinuities are present 
(see e.g. [26| and references therein) . These methods use 
the equations in conservative form together with approx- 
imate or exact Riemann solvers to calculate the numeri- 
cal fluxes between neighboring computational cells. This 
fact guarantees the capture of all discontinuities (e.g., 
shock waves) which naturally appear in spatial solutions 
of the non-linear hyperbolic equations and allows high ac- 
curacy in regions where the fluid flow is smooth. These 
schemes depend on the spectral decomposition of the Ja- 
cobian matrix dT % /dU which has a complete set of lin- 
early independent eigenvectors (characteristic fields) r^, 
and corresponding eigenvalues (characteristic velocities) 
such that 



OF* 



dU 



[r;] = Xi[Ti], 



o,- 



(6) 



In the one dimensional case and in a flat Universe the 
eigenvalues are given by [U HH 



Ao = v x , 



A± = 



Vx ± C s 
1 ± VC S 



(7) 



Xi and Fj are computed by Eqs. (0 — using the follow- 
ing weighted quantities: 



VWri 



i p r v r 



Pi 
plhi 



fPr 



(11) 
(12) 

(13) 



The quantities { AoJfe} represent the jumps of the charac- 
teristic variables across each characteristic field and are 
obtained from the inversion of: 



AU = U r - U, 



3 

E 



(14) 



In order to avoid spurious oscillations we use a mono- 
tone upstream centered scheme for conservation laws 
(MUSCL), with a standard "minmod" sloper limit [32[ 
for the reconstruction of the cell centered quantities be- 
fore the computation of the numerical fluxes. The inte- 
gration in time is performed by using a third order strong- 
stability-preserving Runge-Kutta scheme with five stages 
331. 



IV. GRAVITATIONAL RADIATION SPECTRA 



that correspond to the slopes of three families of curves, 
known as material (Ao) and acoustic (A±) waves. The 
rclativistic speed of sound c s for an EoS in which the 
pressure p is a function of p and e is given by 



where x 



dp 



dp 
~dE 

dp 



x 
h 



p K 



, S is the entropy per particle, 



and E = p(l + e) the total energy density. 

The linearly independent eigenvectors corresponding 
to each eigenvalue are 



ro = 



hW(K - pc? s ) 



V x , 1 



KW{k- pel) 



l,hW\ 



1 



i T 



±- 



-,hW- 



1 



(8) 
(9) 



1 — v\± 1 — v\± 

Our numerical code is based on the Godunov method 
[29| and the numerical fluxes are obtained using the well 
known Riemann solver of Roe 130, 31]: 



fi+1/2 — 



fl)-\ >^ r k .X A„-/, . 1 Id) 



fe=o,± 



which has the advantage that can be straightforwardly 
implemented for an arbitrary EoS. The numerical fluxes 
f) and f r are computed using respectively the primitive 
variables to the left and right of the i + | interface. Both 



To determine the gravitational signal emitted by the 
motion of the primordial fluid, we adopt the formalism 
presented in Ref. (23|, which has also been used in Refs. 
[Til IT2I . among others]. All the information needed to 
calculate the gravitational wave spectrum is contained 
in purely spatial components of the stress-energy tensor 
T )iv (x,t) that, in our case, is obtained from the hydro- 
dynamic simulation. Following this treatment, the grav- 
itational energy radiated per solid angle is given by 

— =2GA ijtlm (k) J uj 2 T^(k,u)T lm (k,u)diu, (15) 

where 

A-ij,lm{k) — fiil^jm 2kjk nl 5n + ^kikjkikm 

is know as projection tensor. The stress-energy tensor 
T^(k,uj) is obtained by a fast Fourier transform (FFT) 
defined as fil [H 



T y (fe,w) = ^- J J d 3 xdt T l3 {x,t)e-^ k - x - ut \ (16) 

Since our hydrodynamic simulations are performed in 
one spatial dimension, our problem is axially symmetric 
about the z axis. Thus, we can take without loss of 
generality [34{ 



sin0, k y — 0, k z 



cos 9. 



(17) 
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Therefore, the projector tensor reads: 



(18a) 
(18b) 

(18c) 



-^2 

A-ijilm 3iz$jz3lz3mz — ^zz,zz = ^(l ~~ ^ z) 



\ sin 4 (9. (18d) 



Finally, the total energy per unit frequency interval 
emitted as gravitational radiation is given by (c.f. [34j ) 



dE 

(lid 



= Glo 2 J |T^(fe,w)sin 2 6» 



+T xx (k,cu) cos 2 9- T yy (k,uj)\ dVt. (19) 



where T zz (k,tu), T xx (k,u) and T»(fe,u) are the Fourier 
transforms of 



T zz (x,t) = hpv 2 z -p, 



T xx {x,t) = T yy {x,t)=p. 



(20) 



(21) 



Integrating Eq. (fT9"|) over the solid angle, the same ex- 
pression as in Eq. (23) of Ref. [n[ is obtained 

g = ^G^ 2 |T"(fc,u,) - T»*(fc, W )| 2 . (22) 

In order to compare our results with eLISA/NGO, we 
describe the spectrum in terms of a characteristic ampli- 
tude of the stochastic background defined as [IBJ 



1Hz 



h c (f) = 1.3 x KT^ficvK/W I — ) , (23) 

where / = w/(2?r), ft = Ho/(100 km s' 1 Mpc -1 ) being 
Hq the Hubble constant, and flow is the energy density 
of the gravitational waves, 



^gw = — (f~u ) J 

Pc V d f 



(24) 



being p c = 3Hq/8itG the critical density. Finally, we 
have to consider the redshift suffered by the waves in 
their way to the present Universe [35| : 



100 



1/3 



Q GW = 1.67 x 10 _5 /io 2 I — J CW*, (26) 

where the subscript corresponds to present values and 
the subscript * to the values at the epoch of the transi- 
tion. 



TURBULENCE IN THE CROSSOVER QCD 
TRANSITION 



Since we are considering a crossover transition at the 
QCD epoch, we don't expect large perturbations near the 
critical temperature as it would be the case for a first or- 
der transition. Thus, we assume that fluctuations present 
at the QCD epoch were generated by some event in the 
previous history of the Universe, e.g. at the electroweak 
phase transition that occurred at t ~ 10 2 s — 10 _10 s. 
Hence, one may consider that the size of the larger fluc- 
tuations at the QCD era are of the order of the size 
of the horizon at the electroweak era; i.e. we have 



clew = O-QCD 



x t 



1/2 
EW 



X t 



-1/2 
QCD 



lm. In other words, 



the injection of energy happens at the electroweak scale 
(~ 1 m) or at an even smaller scale related to inflation 
[5l| . cosmic strings, etc. Q. This fluctuations are con- 
jectured to survive until the beginning of the QCD phase 
transition due to the extremely large Reynolds number 
of the primordial fluid (3(| . 

In turbulent fluids, we usually have a large stirring 
length scale L at which turbulent eddies are injected in 
the fluid. As said above, we are assuming that large 
fluctuations occur at (or before) the electroweak era and 
therefore, the stirring scale is L < lm. As the larger 
eddies break down to smaller ones, there is a cascade 
of kinetic energy from large to small scales. This cas- 
cade effect will stop at a small damping length scale Id 
determined by the viscosity of the fluid. Roughly, the 
damping scale is determined by a Reynolds number close 
to 1, Re ~ 1. Thus we can estimate the order of magni- 
tude of the damping scale as Id ~ Tlev/v, where v is the 
kinematic viscosity of the fluid and v its typical velocity. 

The shear viscosity r\ = v(E + p) of matter at the 
QCD epoch can be derived from heavy-ion collision ex- 
periments at the Relativistic Heavy-Ion Collider (RHIC) 
and the Large Hadron Collider (LHC). In fact, these 
experiments show that the quark-gluon plasma at tem- 
peratures T c < T < 2T C behaves as a nearly perfect 
quark-gluon liquid with a viscosity per entropy den- 
sity in the range 1 < 4n(ri/s)QGP < 2.5 [13, HH, i.e. 
approaching the Kovtun-Policastro-Son-Starinets lower 
bound r\js > l/(47r) (39|. Since the baryon chemical 
potential and temperature of the RHIC/LHC fluid are 
not exceedingly different to what it is expected at the 
early Universe we may assume that the primordial fluid 
at the QCD transition behaves as a perfect fluid as well. 
According to this, we would have Id ~ 0.1 fm which is 14 
orders of magnitude smaller than the size of the cell in our 
numerical simulations, which is Ax — 0.02 m. However, 
the contribution of leptons and photons to the viscosity 
may be significant. In fact, theoretical estimations [40( 
of the viscosity of the quark-gluon plasma (QGP) in the 
early Universe give values that are roughly one order of 
magnitude larger than the experimental values obtained 
at RHIC and LHC. The viscosity may be up to three or- 
ders of magnitude larger if electrons, muons and photons 
are included in the QGP, and even eight orders of mag- 
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nitude larger if neutrinos are considered [40]. In spite 
of this, the viscosity is still very small and the Reynolds 
number very large, as generally considered in the liter- 
ature. Actually, when the most viscous composition is 
considered, the damping length scale is Id < 1CP 7 m, 
which is still orders of magnitude smaller than the size of 
the cell in our numerical simulations. Therefore, the use 
of the ideal hydrodynamic equations presented in Sec. [TT] 
is fully justified. 

Up to now, most works studying the QCD transition 
have considered that the distribution of turbulent kinetic 
energy density has a stationary spectrum consistent with 
the Kolmogorov phenomenology fl2i 13,_41|. i.e. scal- 
ing as ~ fc -2 / 3 or variations of this |l4l. 42j|. This spec- 
trum is very frequent in astrophysical environments, like 
the interplanetary, interstellar or intergalactic medium, 
which have a significant viscosity and a constant injec- 
tion of energy at the large scales (induced e.g. by solar 
flares, supernovae, etc.). This energy is transferred to 
small scales through a cascading effect until ~ l D ; where 
kinetic energy is dissipated into heat due to the fluid vis- 
cosity, resulting a spectrum with a negative slope. In the 
cosmological case, the assumption of a Kolmogorov-like 
spectrum may be acceptable when there is a continuous 
injection of energy, as in the case of a first order phase 
transition, for which bubble collisions introduce a con- 
tinuous stirring. However, this stirring is not present in 
a crossover transition, and therefore the assumption of 
a Kolmogorov spectrum is unjustified. Moreover, since 
the viscosity of the primordial fluid is tiny, dissipation at 
small scales does not occur at the same rate as the energy 
is accumulated by cascading. Therefore, a very different 
turbulent spectrum is to be expected. 



VI. RESULTS 

We have carried out relativistic one-dimensional hydro- 
dynamic simulations for an ideal non-viscous fluid em- 
ploying the fluid equations and EoS of Sec. [TTJ We con- 
sider a computational domain with a length of 100 m, 
with 16.384 spatial cells (2 14 ) and evolve the system for 
times larger than 1 ^s. The choice of these intervals is re- 
lated to the band of the spectrum of gravitational waves 
that we want to obtain from the numerical simulations, 
i.e. ~ 10~ 5 — 10~ 3 Hz according to the eLISA's sensitivity 
curve (see below). We consider reflective boundary con- 
ditions, so that we may see the Universe as a set of var- 
ious contiguous domains of 100 m. Due to choice of the 
boundary conditions, disturbances are reflected several 
times during the simulation, representing the fluid inter- 
action with neighboring regions having similar profiles. 
In agreement with the above discussion, turbulence is in- 
cluded only through inhomogeneities in the initial condi- 
tion. We have considered three different kinds of random 
initial profiles: (a) random temperature inhomogeneities 
in a fluid at rest, (b) random velocity fluctuations within 
a fluid with an initial uniform temperature, and (c) ran- 
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Figure 2: Top: Initial temperature profile with random tem- 
perature fluctuations of maximum amplitude (AT)/T C = 
10~~ 2 . The fluid is considered at rest at t = 0. Center: Tem- 
perature profile after 1/is of evolution. Temperature inho- 
mogeneities are smoothed and the fluid develops a turbulent 
motion. Bottom: Velocity spectrum (v 2 (k)) of the turbulent 
motion of the fluid at t = 1/j,s. There is an energy cascading 
from the larger to the smaller scales. Energy is not dissipated 
at the smallest scale because of the negligible viscosity of the 
primordial fluid. We have considered several different random 
initial conditions for the temperature profile, all with maxi- 
mum amplitude (AT)/T C = 10~ 2 . All of them present the 
same behavior as presented in this figure. 



dom temperature and velocity fluctuations. Since the 
size of possible disturbances at the time of the transi- 
tion is unknown we consider fluctuations of maximum 
amplitude AT/T C ~ 10~ 2 , 10~ 3 , 10~ 4 around T c = 170 
MeV and/or Av/c ~ 10 _1 , 10~ 2 , 10 -3 in the initial pro- 
file. Our goal is to determine the smallest fluctuation 
amplitude that would be detected by eLISA, considering 
the motion of the fluid induced by the initial condition 
as a source of gravitational radiation. The hydrodynamic 
simulations provide the tensor T(x, t) from which we ob- 
tain T(k,u>) through a FFT. This allows the calculation 



6 




/[Hz] 

Figure 3: Spectra of gravitational waves for hydrodynamic 
simulations with different initial conditions at the beginning 
of the QCD epoch. We considered random temperature fluc- 
tuations of maximum amplitude (AT)/T C = 1CT 2 (see Fig. [2J 
as well as 10 -3 and 10 -4 . For comparison we show the sen- 
sitivity curve of eLISA/NGO computed using the expected 
instrumental noise and the confusion noise generated by un- 
resolved galactic binaries [221 ]. 



of the spectrum of the gravitational radiation emitted by 
the fluid through Eq. (|2"2"j) . 

In Fig. [2] we consider random temperature inhomo- 
geneities with a maximum amplitude of 10 -2 in a fluid 
that is initially at rest. The initial temperature profile 
(top panel) presents peaks and valleys of different sizes. 
As the fluid evolves, these temperature gradients induce 
the motion of the fluid, temperature inhomogeneities are 
smoothed (see central panel of Fig. [3]) , and the fluid de- 
velops a turbulent motion with velocities up to ~ 0.02c. 
As discussed before, there is a cascading that transports 
kinetic energy to the smallest scales where there is no dis- 
sipation due to the absence of viscosity. This behavior is 
apparent in the bottom panel of Fig. [21 which shows that 
a large amount of kinetic energy is accumulated at the 
smallest scales (large k) at the end of the simulation. We 
performed similar simulations with random temperature 
inhomogeneities with a maximum amplitude of 10 -3 and 
10~ 4 . The hydrodynamic evolution is not shown here 
because the behavior is qualitatively the same. 

In Fig. [3] we can see a comparison between the gravi- 
tational wave spectra arising from the three initial tem- 
perature gradients considered here. In the case with 
AT/T C < 10~ 2 the signal would be detected by eLISA 
for a wide range of frequencies. For the simulation with 
AT/T C < 10~ 3 the signal is entirely below but not too 
far from the eLISA's threshold. In fact, we can show that 
the signal for fluctuations of about AT/T C <3x 10~ 3 is 
partially overlapping the eLISA's sensitivity curve. Fluc- 
tuations with AT/T C < 10 don't lead to enough motion 
in the fluid to emit a significant amount of gravitational 
radiation. 

We also considered an initial condition with constant 
temperature T — 170 MeV and a random velocity distri- 
bution. For an initial maximum amplitude of (Av)/c = 
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Figure 4: Top: Initial velocity profile with random veloc- 
ity fluctuations of maximum amplitude (Av)/c = 10~ 2 and 
constant temperature T = T c = 170 MeV. Center: Velocity 
profile after 1 /is of evolution. Bottom: Velocity spectrum 
(v 2 (k)} at t = and t = lfxs. The final spectrum is steeper 
due to the energy cascading from the large to the small scales. 
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Figure 5: Spectrum of gravitational waves for simulations 
with initial random velocities (Au)/c = 10 -3 , 10 -2 , 10 _1 . 
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Figure 6: Left panels: Initial profiles with random velocity and temperature fluctuations of maximum amplitude (Av)/c w 10 -2 
and (AT)/T C = 10 -2 . Central panels: velocity and temperature profiles after 1/is of evolution. Right panels: same as before 
but after 10/xs of evolution. 
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Figure 7: Velocity spectrum at different times for the initial 
conditions of the previous figure. 
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Figure 8: Spectrum of gravitational waves for (Av)/c = 10 
and different initial temperature fluctuations. 



0.01 (see top panel of Fig. @| the turbulent motion creates 
temperature gradients of the order of AT/T C w 3 x 10~ 3 
which at the end of simulation fall to about AT/T C « 
1.5 x 10 -3 with maximum velocities roughly the half of 
the initial (see central panel of Fig. 2]). On the bottom 
panel of Fig. HJ we show the initial and final velocity spec- 
tra, which present the same behavior as in the previous 
simulations, confirming the expectation that the spectra 
must be very different to the Kolmogorov power law. For 
maximum initial velocities of (Av)/c ^0.1, the induced 
temperature gradients are an order of magnitude larger 
than in the previous case and the strong motion of the 
fluid leads to a large amount of gravitational radiation 
(Fig. El). Notice that the spectrum with (Av)/c < 10~ 2 
is almost overlapping the eLISA sensitivity curve for fre- 
quencies larger that 10 -4 Hz as in the simulation with 
initial condition AT/T C < 3 x 10~ 2 , (Av)/c = 0. In spite 
of these initial conditions being very different, the system 
evolves in both cases to similar temperature and veloc- 
ity profiles in a timescale considerably smaller than the 
total time of the simulation. In fact, simulations with 
concomitant velocity and temperature gradients in the 
initial condition give essentially the same results as be- 
fore, i.e. a progressive homogenization of both the tem- 
perature and velocity profiles, an energy cascade to the 
smaller scales, and a gravitational wave spectrum with 
a peak around 10 -4 Hz (see Figs. |51-[5]). This peak in 
the spectrum is due to an exponential decay of the fluid 
velocity. Actually, we verified that \v\ ~ exp(— t/r) with 
t ~ 10~ 8 s for most values of the wave number k, showing 
that most of the fluid motion happens within this time 
interval. Since gravitational wave emission is associated 
with the level of turbulent motion in the fluid, we expect 



8 



a maximum in the spectrum at a frequency ~ r _1 , which 
from Eq. (|25D turns out to be ~ 10~ 4 Hz; i.e. consistent 
with the maximum in the spectra of Figs. [3j [5] and [8] 

VII. SUMMARY AND CONCLUSIONS 

In this paper, we have studied the evolution of primor- 
dial turbulence in the universe at the QCD epoch using 
a state-of-the-art equation of state derived from lattice 
QCD simulations [3[. Since the transition is a crossover, 
we don't expect large perturbations near the critical tem- 
perature as it would be the case for a first order transi- 
tion. Thus, we assume that temperature and velocity 
fluctuations were generated by some event in the previ- 
ous history of the Universe (e.g. a first order transition 
at the electroweak scale or at an even smaller scale re- 
lated to inflation, cosmic strings, etc.) with typical sizes 
smaller than the size of the horizon at the electroweak 
era, L < lm. Due to the extremely large Reynolds num- 
ber of the primordial fluid we consider that these inho- 
mogeneities are able to survive until the QCD epoch. 

The primordial fluid at the QCD epoch is assumed to 
be non-viscous, based on the fact that the viscosity per 
entropy density of the quark gluon plasma obtained from 
heavy-ion collision experiments at the RHIC and the 
LHC is extremely small, more specifically, in the range 
1 < A-k(t]/s)qgp < 2.5 at temperatures T c < T < 2T C 
(37l [Hj]. For the hydrodynamic simulations, we con- 



sidered a one dimensional computational domain with 
a length of 100 m divided in 2 14 cells, injected random 
temperature and velocity fluctuations in the initial condi- 
tions, and evolved the system for times larger than 1 /is. 
Our results show that the velocity spectrum is very differ- 
ent from the Kolmogorov power law considered in most 
studies of primordial turbulence that focus on first order 
transitions. This is due to the fact that there is no contin- 
uous injection of energy in the system and the viscosity 
of the fluid is negligible. Thus, as kinetic energy cascades 
from the larger to the smaller scales, a large amount of 
kinetic energy is accumulated at the smallest scales due 
to the lack of dissipation (see bottom panels of Figs. [5] 
andU and Fig. E}. 

We have obtained the spectrum of the gravitational 
radiation emitted by the motion of the fluid for differ- 
ent initial profiles that include random temperature and 
velocity fluctuations of different maximum amplitudes. 
We find that if typical fluctuations have an amplitude 
{Av)/c > 10~ 2 and/or AT/T C > 10~ 3 , they would be 
detected by eLISA at frequencies larger than ~ 10 ~ 4 Hz. 
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